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Abstract 


In this paper we present a one-dimensional transient model for the membrane electrode assembly of a polymer-electrolyte fuel cell. In earlier work 
we established a framework to describe the water balance in a steady-state, non-isothermal cathode model that explicitly included an agglomerate 
catalyst layer component. This paper extends that work in several directions, explicitly incorporating components of the anode, including a micro- 
porous layer, and accounting for electronic potential variations, gas convection and time dependance. The inclusion of temperature effects, which 
are vital to the correct description of condensation and evaporation, is new to transient modelling. Several examples of the modelling results are 
given in the form of potentiostatic sweeps and compared to experimental results. Excellent qualitative agreement is demonstrated, particularly in 
regard to the phenomenon of hysteresis, a manifestation of the sensitive response of the system to the presence of water. Results pertaining to pore 
size, contact angle and the presence of a micro-porous layer are presented and future work is discussed. 


© 2006 Published by Elsevier B.V. 
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1. Introduction 


The proton exchange membrane fuel cell (PEMFC) is a 
promising source of clean and efficient energy, particularly 
for application in the automotive industry. As with other fuel 
cells, the basic principle of the PEMFC is to convert chem- 
ical into electrical energy. This aim is achieved by oxidation 
of hydrogen at the anode and simultaneous oxygen reduction 
at the cathode, with a net production of water. The reactants 
are delivered through a graphite paper, termed the gas dif- 
fusion layer (GDL), and reaction takes place in thin regions 
termed the catalyst layers. The catalyst in question is plat- 
inum (Pt), which is required to enhance the rate of reaction 
by providing an alternative pathway to the final products. 
Clearly, optimal use of the Pt is a key concern in all PEMFC 
designs. 

The catalyst layers are perhaps the most complex compo- 
nents of the fuel cell. Through the solid components (carbon 
grains, supporting smaller catalyst particles, and electrolyte), 
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flow water vapour and components of air. Liquid water resides 
in the gas pores and absorbed water within the electrolyte. The 
electrons migrate through the carbon components of the layer 
while the protons migrate through pathways provided by the 
electrolyte—typically Nafion. The hydrophilic regions in the 
cathode catalyst layer (CCL) can cause retention of the water 
produced during reaction, which restricts the ingress of oxy- 
gen. Moreover, reaction is limited by the availability of Pt, and 
can only occur at points of contact between the Pt, carbon and 
electrolyte. 

The properties of the membrane and the GDL, with possi- 
bly a micro-porous layer, also play a central role in the control 
of water. It is essential that the protons generated at the anode 
reach the cathode freely, requiring that the membrane remain 
well hydrated. The GDL must allow sufficient oxygen from 
the cathode channel to enter the CCL, but must be sufficiently 
hydrophobic to expel any build-up of liquid water beyond that 
required for optimal hydration of the membrane. The importance 
of water management explains the emphasis placed on capturing 
the effects of (particularly) liquid water in modelling studies [1— 
9]. However, a lack of understanding of the interaction between, 
and transport of, the three phases of water continues to retard 
progress towards a predictive model. 
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Steady-state PEMFC models are numerous. Because of the 
many challenges faced with the task of modelling the com- 
plex physics and geometry, a host of simplifying assumptions 
are regularly adopted, typically involving the catalyst layers, 
liquid water, electron transport, thermal energy and time depen- 
dence. For example, Natarajan and Nguyen [5], use a simplified 
and regularized form of the liquid water permeability in their 
isothermal model and do not explicitly model the membrane, 
catalyst layer and anode. Siegel et al. [6], incorporate the CCL 
and membrane, making the assumption that water exists only as 
vapour or as a species dissolved in the electrolyte/‘membrane, 
and are not able to account for liquid water effects directly. 
The same is true of the model in [10], which treats all three 
forms of water as a single phase. In [7], all three forms are 
included but with a highly regularized capillary diffusion coef- 
ficient. In [11-14], and references therein, the authors describe 
the water balance using the so-called multi-phase mixture (M7) 
model, developed in [15], where thermodynamic equilibrium 
between liquid water and vapour is assumed to prevail through- 
out (saturated gas phase and no vapour diffusion) and saturation 
is calculated a-posteriori based on the equilibrium assumption 
and mixture water concentration. The advantage of this approach 
lies in its simplification of the multi-phase problem. Models 
employing the M? approach have until recently assumed isother- 
mal conditions, [12]. In this paper the M? simplification is not 
used. 

In contrast, very few transient models of PEMFC exist, 
despite the obvious transient nature of the cell operation in 
the primary application of powering automobiles. Um et al. 
developed a single-phase, isothermal 2D transient model in 
[16] and a similar 3D model in [17], though transient simu- 
lation results are not presented in the latter. 3D models can 
also be found in [18-21], none of which include liquid water 
transport and an energy balance—vital to describing flooding, 
membrane hydration and phase change. A more comprehen- 
sive one-dimensional model is presented in [22] in which liquid 
water is accounted for in a fully two-phase manner. How- 
ever, the model is isothermal and in comparisons with potential 
sweep experiments it fails to capture the sensitive balance 
between flooding effects and membrane hydration. For a more 
detailed review of fuel cell modelling to 2004 we refer to 
[23]. 

The aim of this paper is to reproduce the phenomenon 
observed in the aforementioned experiments and those in [24] 
by developing a transient, non-isothermal model that explic- 
itly incorporates the GDL, catalyst layers and membrane. In the 
present case, a one-dimensional formulation, based largely on 
our earlier work in [25], is shown to be sufficient to qualita- 
tively capture observed behaviour — by comparison with [22,24] 
— and predict performance. Efforts to extend the model to two 
and three dimensions are ongoing. 


2. Model assumptions and equations 


The various features and assumptions that form the basis of 
our model are now listed. 
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Fig. 1. Model geometry. 


1. One-dimensional full MEA transient. The domain includes 
the entire membrane electrode assembly from the interface 
between the gas channels and GDL to the membrane (see 
Fig. 1). Each component is modelled explicitly. We allow for 
the inclusion of micro-porous layers between the GDL and 
catalyst layers. These layers serve two primary purposes: to 
prevent flooding in the catalyst layers and to decrease the 
contact resistance between the catalyst layer and diffusion 
medium [26,27]. 

2. Catalyst layers. We assume that the carbon support forms 
spherical clusters (agglomerates). Surrounding each agglom- 
erate is an electrolyte layer, on which a layer of liquid water 
can also exist. We assume that the contact between the carbon 
agglomerates is sufficient to ensure the free flow of protons 
and electrons, and that the connecting paths posses a neg- 
ligible volume fraction (a Bruggemann correction is used 
to account for the tortuosity of the flow paths). Within the 
agglomerates exists a fraction of the electrolyte, which we 
assume forms a negligible surface area of contact with any 
Pt that may also exist inside the agglomerates. It is possible 
to derive an effectiveness factor to remove this assumption. 
However, since it is unclear that significant reaction occurs 
within the agglomerates [28-31], an effectiveness factor is 
not employed. The pores between agglomerates are referred 
to as primary pores, distinct from the smaller pores between 
the carbon particles. 

3. Reactant transport. Diffusion in the primary pores of the 
catalyst layers is assumed to be predominantly continuum, 
given their characteristic size and the typical operating pres- 
sures. Assuming that the oxygen and water vapour are small 
concentrations in a nitrogen gas, we employ Fick’s law on 
the cathode side. On the anode side we assume the typi- 
cal composition of hydrogen and vapour, with hydrogen the 
dominating component. The oxygen and hydrogen in the pri- 
mary pores of the catalyst layers reach the Pt particles on the 
agglomerate surfaces by dissolving in the water/electrolyte 
and diffusing across the surrounding layers. Henry’s law is 
used to describe the equilibrium in oxygen concentration at 
the interface between the gas and liquid/solid phases and 
mass is allowed to be exchanged freely at the surfaces. That 
is, the deviation from this equilibrium provides a driving 
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Table 1 
Sources and sinks for the vapour equation (3) 
ACL CCL GDL/MPL Meaning 
See — heel Kap — Psat) —Npc(X yp — Psat) —hpe(Xvp — Psat) Condensation/evaporation 
Sad hay(Ca — C3) hav(Ca — C3) 0 H20"4/H, 0's exchange 


force for interfacial mass transfer, which occurs at a finite 
rate. 

4. Proton and electron concentration. The protons are assumed 
to be transported in the form of hydronium ions, H3O*. 
We explicitly model proton and electron transport through 
the ionomer, membrane, carbon support and GDL, assum- 
ing electro-neutrality. As in [16,32,22,20,21], we assume 
a pseudo steady state for the electron and proton trans- 
port/generation. The justification for this can be found in 
[20]. 

5. Water. The model accounts for water in all three forms; as 
a dissolved species, as vapour and as liquid. Liquid water 
resides in the primary pores. We assume that the net water 
produced is in dissolved form, consistent with the reac- 
tion location and the strongly hydrophillic nature of the 
ionomer. Condensation and evaporation are modelled using 
the approach in [2,4,5], and references therein. In this rep- 
resentation the transfer occurs at a finite rate, dictated by 
the deviation of the local thermodynamic state from equi- 
librium. In a similar fashion, we introduce phase change 
between vapour and dissolved water and between liquid and 
dissolved water by considering the deviation from an equi- 
librium between the phases. Details are presented later. 

6. Convective flow. As in [5,22,33] the steady-state Darcy’s law 
is employed for the convective flow of gas and liquid in the 
pores. This approach is also found in mining, hydrology and 
combustion literature [34], where a steady-state Darcy’s law 
is a standard approximation. 

7. Temperature. The model is non-isothermal, which is neces- 
sary to model phase change accurately and to study the effects 
of water activity. We do however assume a single temperature 
for all phases. 


2.1. Gas phase equations 


Let s denote liquid water (H2 old) saturation, Co, Cn, Cy and 
Cy, respectively the concentrations of oxygen, nitrogen, hydro- 
gen and water vapour in the pore space, and Ce, i € {O, H}, 
the concentrations of oxygen and hydrogen in the catalyst layer 
ionomer. The subscript j refers to anode (a) and cathode (c), and 
the subscript i refers to oxygen (O) and hydrogen (H). Taking 
into account reaction, diffusion and absorption, mass balances 


yield 
d d aC; 
J rat dr (o, F wei 
= —hpe,i(HiCi — Cie) d) 
d d dCN 
1 D 5 = 2 
ezk 5s)CN] Ox ( N= EN 0 (2) 


d d d, 
€ dr DI = sev] ox G Ox 
where Dp,i, DN and Dy, i € {O,H}, are the diffusion coeffi- 
cients of oxygen, hydrogen, nitrogen and water vapour (H2O0%"?) 
through the GDL and primary pore of the catalyst layers. They 
are subject to a Bruggeman correction and their dependence on 
temperature and pressure follows the Chapman—Enskog formula 
[35], shown in Table 5. In these expressions T is temperature and 
P is the gas pressure. The porosity e takes the value e in the 
catalyst layers, €g in the GDL and €mp in the micro-porous layers. 

The quantity vg is the gas velocity, assumed to follow Darcy’s 
law 


H = See + VSad (3) 


K oP 
v = —— 01 sy —, 
g u ( s) ax 
where v is the permeability of the GDL, micro-porous or catalyst 
layers and u is the dynamic viscosity of the gas. The permeability 


can be approximated by the Kozeny—Carman law, 


(4) 


Ê ei 


= E EEN (5) 


where d is a mean pore diameter and K is the so-called Kozeny— 
Carman constant. The two source terms in Eq. (3) are defined in 
Table 1. 

The quantity v = 1800 mol m7? is the fixed-charge site con- 
centration of the membrane; Ru: T, C Gel are the reaction rates 
(based on the Butler—Volmer law); Ree: is a volumetric mass 
transfer coefficient from gas to electrolyte (on the air side) and 
H; is a dimensionless Henry constant. The reaction kinetics at 
the two electrodes can be summarized as follows: 


HOR : Hə — 2H* =2e , 
ORR : 2H)O — O; AHT = Ae (6) 


We denote by vj, j € {a,c}, the stoichiometric coefficients and 
by n the number of electrons transferred. 

The diffusive flux of reactants from the pore space to the 
interface between the gas and electrolyte/water is balanced by 
the amount absorbed into the electrolyte/water. To prescribe the 
mass-transfer coefficient, Mpe,;, we use the formula for the Sher- 
wood number given in [35] for flow past a spherical particle 


Sh =2+0.6Re!/? Sh, (7) 


where Re is the Reynolds number and Sc is the Schmidt num- 
ber. For a predominantly diffusive process, we can approximate 
Sh = 2. From the formula hpe,i/Sa = ShDp,i/D, where D is an 
average pore diameter and sa is the specific surface area of the 
agglomerates, we obtain bn, 7 = O(10°) or greater. We point 
out that estimating mass transfer coefficients is far from a trivial 
exercise and a detailed discussion will not be attempted here. 
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Table 2 
Sources and sinks for the dissolved reactants, potential and dissolved-water equations (10), (12) and (14) 
ACL CCL Membrane GDL/MPL Meaning 
Sex Dee pt Dun — Cy,e) hpe,o(HoCo — Co,e) 0 - Reactant absorption 
Sr —VaRa(na, T, Cy e)/” —VeRe(Ne, T, Co,)/n 0 - Reactant depletion 
So —nFS,/Vq —nFS;/V¢ 0 H30* source/sink 
Sy — So — So - e~ source/sink 
Sw 0 Rate T, Co,.)/(2v) 0 — Water production 
Sat Iya(Ca — C3) Iya(Ca — C3) - H20"9/H20%s exchange 


2.2. Membrane and carbon phases 


The electrolyte volume fraction, €e, will increase as the elec- 
trolyte swells with water. We separate the electrolyte film that 
coats the agglomerates (volume fraction e) from that contained 
in the interiors (volume fraction el) and write e? = ef + ei. The 
Pt inside the agglomerates is assumed to be inactive, that is, there 
is a negligible surface area of Pt in contact with the fraction eL of 
electrolyte. The combined volume fraction of the carbon, Pt and 
small pores, €a, is assumed constant. The volume fraction of pri- 
mary pores is then given by €p = 1 — €a — €e. The swelling of 
the electrolyte is assumed to impact on the thickness of the elec- 
trolyte film. To incorporate it we use the following approximate 
relationship 


ee = el + ksà, (8) 


where ks = 0.0126, and A is the membrane water content (mol 
H20/mol SO; 1. given by 


Ca 


Ass 9 
1 — k;Ca ©) 


where Ca is the dissolved water (H208) concentration in the 
electrolyte, normalized with respect to (e, divided by) v. 

The equations for reactant concentration in the ionomer and 
membrane are given by 


OC; d OC; 
oe Lon = Sr + Sex, (10) 
where e = ee in the catalyst layers and e = 1 in the membrane. 
Note that we use a Bruggeman correction. The diffusion coeffi- 
cient of oxygen through the electrolyte has the following form 
(taken form [36]) 


Deo = 3.1 x 10, (11) 


For hydrogen we use the constant value in Table 5. 

To derive equations for the potentials in the membrane 
and carbon phases, @ and y respectively, we consider a mass 
balance for H3O* and electrons, with the assumptions of electro- 
neutrality and steady-state 


ox (< a ll g dr (< Ge ox eae 


(12) 


where de and o; are the protonic and (effective) electronic con- 
ductivity respectively and F is Faraday’s constant. Note that we 
have used a Bruggeman correction in both equations. € = €e 
in the catalyst layers and e = 1 in the membrane for the ionic 


potential and € = 1 — €g (€ = 1 — €mp) in the GDL (MPL) and 
€ = €c in the catalyst layers for the electronic potential. The 
source terms are given in Table 2. 

For ce we use the law derived in [37], 


1 1 
Oe = exp (1286 E di ell (0.5144 — 0.326). (13) 


The algebraic term on the right-hand side of (13) accounts for the 
effect of hydration on the conductivity. It is based on empirical 
data collected by Springer et al. [37]. Values for os have been 
discussed by Sun and Karan in [38]. To be consistent with their 
interpretation, we use the value o, = 500S ml, 

The mass balance for water dissolved in the electrolyte is 
written as follows 


aa A (py Ca, 2 a) _ 
Ca ae ge aa a) 


Sad + Sw — Sai, 


(14) 


in which e = 1 for the membrane and e = ee for the catalyst 
layers. The term Sw represents water production and Sq repre- 
sents mass transfer between the liquid and dissolved phases (see 
Section 2.5). Dg is the diffusion coefficient for dissolved water 
[39] 


Dg = 3.1 x 1077 A (e284 — 1) e72436/T 


Da = 4.17 x IO BA + 161 e~*) e 2436/7 


(0 <A <3), 
(3 <à <22). 
(15) 


The source term Saq is defined in Table 1 


2.3. Energy 


An equation for the temperature, T, is derived from an energy 
balance of conduction, convective heat flux and heat sources 


a 
Ge (espCiT + €(1 — s)pgCgT + psCsT) 


ox 
= Sact + Srev + Sohm + Spc, (16) 


d d oT 
+ ER Lean? + eil — S)pgCgT) — — EA 


where o (pg) is the density of liquid water (the gas phase), Cı 
(Cg) is the specific heat capacity of liquid water (the gas phase), 
and k is the effective thermal conductivity, found from a volume 
average of the individual conductivities 


k = kp = S)E + Ke€e + zer + Ka€a, (17) 
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Table 3 
Sources and sinks for the energy equation (16) 

Membrane Catalyst layers GDL/MPL Meaning 
Sact 0 FnjRj(nj, T, Che) 0 Activation losses 
Srey 0 —AjTRi(nj, T, Cf )/n 0 Heat of reaction 
Soim oe(Ob/ Ax)? l olap + eil oya (1 — €g /€mp)’ 20s (04 /3x)? Ohmic losses 
Spe 0 hgihpe(Xvp — Psat) hgih sc Kap — Psat) Heat of evaporation 


where kp, ke, ka and k; are the thermal conductivities of the pore 
space, electrolyte, agglomerates (averaged) and liquid water 
respectively. In the GDL e€ = ég, in the micro-porous layers 
e = €mp and in the catalyst layers e = ép. 

The quantity u is the liquid water interstitial velocity, which 
we define later. The source terms in equation (16) are defined in 
Table 3. In these expressions — A ; is the entropy associated with 
ORR (j = c) and HOR (j = a), hpc is amass transfer coefficient 
(defined below) and hg is the liquid—gas enthalpy change for 
water. 


2.4. Liquid water 


For liquid water, the mass balance is 


eo OS d (ee 
W, ot dr Wi 


n) = Sce + vSal, (18) 


where the source terms, defined in Tables 1 and 2, represent 
condensation/evaporation and interfacial mass transfer between 
the liquid and dissolved phases (the last term on the right-hand 
side is explained and discussed in Section 2.5). W, is the molar 
mass of liquid water, € = €g in the GDL, e = emp in the micro- 
porous layers and e = en in the catalyst layers. vj is determined 
by the Darcy-law approximation to momentum conservation 


K] dn 
u= — 


e l 19 
Uu Ox GC 


where vi, yı and p are respectively the relative permeability, 
viscosity and pressure of the liquid. 

A great deal of debate surrounds both the form and magni- 
tude of the permeability and capillary pressure. Mazumder [4], 
compares the Leverette model employed in [9] with the empiri- 
cal relationship in [5], in which both of these quantities assume 
entirely different forms, yielding markedly different predicted 
levels of saturation. There are also other models of liquid trans- 
port in porous media, of varying degrees of sophistication. For 
example, Nam and Kaviany [40], scale saturation against an 
immobile value (below which the liquid water is discontinuous 
and therefore does not move), and place the scaled saturation 
inside the Leverette function. All of the forms mentioned can 
be easily implemented in our model but, in the absence of con- 
clusive experimental data, we choose to follow the most widely 
employed form, found in [9]. 

The liquid pressure is given as the difference between the gas 
and capillary pressures 


p=P-— pc, (20) 


so that differentiating Eq. (18) and using Eq. (19) yields 


eo OS o d dp. de OP 
ee ape ee eg 
Wi at pW) ox (<a) | bk n caor Hed 


(21) 


In the catalyst layers the capillary pressure and permeability 
are defined as 


Pe = Ks), KICS) = Kes, (22) 


where «e is the absolute permeability of the catalyst layers and 
Ks) is the Leverette function. The quantity oc is defined as 
follows 


o, = o' cos 6 Se (23) 
Cc 


where 62 is the contact angle for the catalyst layers (0 < 8l < 
7/2) and ol is the surface tension. The grouping 24/€p/Kc is 
then equal to the characteristic capillary radius. 

The Leverette function for a hydrophilic medium, for which 
the wetting phase is the liquid, is given by 


Ks) = 1.417(1 — s) — 2.121 e + 1.26201 en (24) 


In the GDL and micro-porous layers 


E 
Pe = Ogs), K\(s) = Kgs”, og =0' cos 68 SCH 
K 
8 
€ 
Pe = omp s),  K(s) = Kaps a Omp = o cos oP — 
(25) 


where Kg (Kmp) is the absolute permeability of the GDL (micro- 
porous layer) and 68 (02) is the contact angle (7/2 < OÈ < r). 
The Leverette function for a hydrophobic medium, for which 
the wetting phase is the gas phase, is given by 


Ks) = 1.4175 — 2.12s? + 1.2625°. (26) 


The liquid—vapour phase-change term, derived from that 
given in [4] and references therein, is driven by the deviation 
from equilibrium, RTC, — Pyar, where Psat is the saturation pres- 
sure of water and the first term is equivalent to the partial pressure 
of the vapour. The latter is defined by Xy P, where P is the pres- 
sure of the gas phase and Xy is the molar fraction of vapour. 
Using the ideal gas law, P = CRT, Xy is given by 


l (27) 


where C = Cy + Cy + Ci is the molar density of the gas (i = O 
at the cathode and i = H in the anode). A relationship for the 
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saturation pressure is provided in the discussion of the bound- 
ary conditions. In order to distinguish between evaporation and 
condensation, the mass-transfer coefficient hp, depends on the 
sign of the driving force Xyp — Psat and its precise form is 


ge ep(1 — s)Xy 


RT Xyp > Psat 


f (28) 
Xyp < Psat 


where kc and ke are the condensation and evaporation rate con- 
stants, whose values are taken from [2]. We use a continuously- 
differentiable implementation of the following form 


8 keep — ELE, (1 H |Xyp — Si 


d 2RT Xyp — Prat 
4 Ke€pSpl (: |Xyp — rt) ; (29) 
EIN Xyp = Psat 


In a similar fashion, the vapour-dissolved phase-change term 
in Eqs. (3) and (14) is driven by the deviation from equilibrium 
between the vapour and dissolved water, Ca — Ci, where C3 is 
the dissolved concentration at equilibrium. The latter is derived 
from the following law for the equilibrium water content mea- 
sured at 80°C (found in [41]) 


Pai 
~ 1+0.0126.*" 
(30) 
where aw is the water vapour activity. To distinguish between 
desorption and water-uptake of the electrolyte (adsorption), the 
mass-transfer coefficient Aay depends on the sign of the driving 
force Ca — Cj. From the experimental and numerical results 


in [42], we can approximate the coefficients of absorption and 
desorption as follows 


Kl a Cu 
Kal — s)A Ca- Ci > 0 , 


A7 = 0.34 10.8aw — 1602+ 14.103, Ci 


hay = (31) 
where «q and ve are given in Table 4. Their weak dependence on 
temperature variations is neglected. The factor 1 — s accounts 
for the blockage of the pores by the liquid water. To implement 
this approach we use a continuously differentiable form of 


1 Cq — C* 
hay = ni H (1 + SEH 


Ca- E 
o Q—sa(1 ICa = Gal (32) 
K. S 
2” Ca- JI 


2.5. Liquid-dissolved water mass transfer and Schroeder 
paradox 


Itis well known that the water content A depends on the water 
activity, dy = Xv P/ Psat. The precise relationship was correlated 
by Hinatsu et al. in [41], 


à = Ay = 0.34 10.8ay — 16a2, + 14.103, (at353K), (33) 


which has a maximum of A, = 9.1 at equilibrium with vapour. 
However, when the electrolyte is submerged in liquid water its 


equilibrium water content appears to jump discontinuously to a 
higher value, given in the following 


22 


EH 34 
1+ 22k Ge 


A = 22 orequivalently Cy=Cj)= 
as dictated by Eq. (9). In an attempt to capture this anomaly, 
known as Schroeder’s paradox, we introduce the mass-transfer 
term Sq in Eqs. (14) and (18), defined in Table 2. This term is 
decomposed into separate terms for absorption and desorption of 
liquid water to and from the ionomer in the catalyst layers. When 
the liquid-equilibrated water content value Cj | is reached, it is 
assumed that desorption of water from the ionomer as a liquid 
takes place, the magnitude of which is driven by the deviation of 
Ca from CH. The desorption part of the term in Table 2 therefore 


takes the form 
Ca — Ci 
zeu) : (35) 


2 Ca — Cal 


des = =M (: + 
where kaes is the coefficient of desorption and could itself depend 
on saturation and water content but for simplicity is assumed to 
be constant. 

In addition we introduce an absorption term that represents, 
in its presence, the uptake of liquid water by the ionomer, the 
magnitude of which is again assumed to be proportional to the 
deviation of Ca from Car However, it is further assumed that 
the liquid water droplets are not contiguous until the immobile 
saturation, Sx, is reached (by definition of the immobile satura- 
tion). By this approximation we are equivalently assuming that 
at s = s, the agglomerates are entirely coated with liquid water 
and equilibrium can be achieved. We use the value of są = 0.1 
given in [40], and refer to that paper for references related to its 
measurement. The absorption term therefore takes the form 


1 = l Co CH 
E E EE EA E 
2 S — Sy 2 


Ca — Cay 
where kaas is the coefficient of absorption, which, for simplicity, 
is assumed to be constant. Notice that the last term in brackets 
in (36) is zero when Cg > C3 so that it does not contribute to 
desorption. The mass transfer term in Table 2 is therefore defined 


as follows 
k s-—S ICa— Cd 
ha — Kes (14 | «| 1 du 
4 A Sy Ca — CA 
k des [Ca — CG 
+ 1+ ———— 37 
2 ( Ca — Cir G7) 


Again this is implemented in a continuously differentiable form 
to prevent the occurrence of singularities in the numerical com- 
putations. 

The values of Kads and kqes are chosen large enough that if 
Ca approaches Cj ; its value does not overshoot it significantly 
(as with the treatment of condensation and evaporation). These 
values are given in Table 5. A list of other parameters and their 
values can also be found in Tables 4 and 5. 
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Table 4 
Parameter values used in the base case calculations 
Symbol Quantity Size References 
Structural 
L Catalyst layer thickness 25 um 
La Membrane thickness 50 um 
Linp GDL + MPL total thickness 200 um 
Lo MPL thickness 30 um 
D Electrolyte volume fraction in agglomerates 0.2 [30 
Ea Volume fraction of agglomerates 0.4215 [30 
Ec Volume fraction of carbon in catalyst layers 0.2 
Eg Porosity of the GDL 0.74 [47 
Emp Porosity of micro-porous layer 0.3 [47 
Rag Agglomerate radius 0.2 um [30 
55 Electrolyte film thickness w/o swelling 15nm [48 
N Agglomerate density 1.257 x 10!° m~? est. 
Apt Specific surface area of Pt 1000 cm? (mg Pt)! [49 
Mpt Pt loading 0.4 (mg Pt) cm~? [3] 
oi (62) Catalyst-layer (GDL) contact angle 90° (120°) [50 
dg, dun, de GDL, MPL, catalyst-layer pore size 10 um, | um, 2 ym 
Electrochemical 
bref i Cathode (anode) exchange current density 1073 (103) A m~? [51 
ef ? Reference O2 (H2) concentration 0.005 (0.02) mol m~? B1 
Qe Cathodic transfer coefficient 0.55 
Qa Anodic transfer coefficient 0.45 
Eo Open circuit potential 0.95 V 
— Ae Entropy associated with ORR 163.7J mol“! K7! [52] 
— ^a Entropy associated with HOR? 0Jmol~! K7! [52] 
Channel conditions 
Te (Ta) Cathode (anode) channel temperature 60°C (60°C) 
dye (awa) Cathode (anode) channel water activity 0.8 (0.8) 
Ce (Ca) Cathode (anode) channel gas concentration P./RTz (Pa/RT,) mol m7 
Cox Oxygen concentration in cathode channel 0.21(C, — Cy.) mol m~? 
CHa Hydrogen concentration in anode channel (Ča — Cy,,) mol m~? 
Po (Pa) Gas pressure in the cathode (anode) channel 30 psig (206.842 kPa) 
Cie Vapour concentration in cathode channel RT PsatrcCe/ Pe mol m~? 
Cra Vapour concentration in anode channel RT, Poata Ca / Pa mol m`? 
Qj Liquid water removal constants 0.075 m7! 


2.6. Reaction rate and limiting current density 


The reaction rates (in mol m~? s7!) are given by the Butler- 
Volmer law and first-order kinetics in reactant concentration 


Kan, T, Cie) = 


Set Ze Ges _ eee) 3 


(38) 


where i = O and j = c ori = Hand j = a. The quantities iref,i 


are the cathode and anode exchange current densities, œa and 
Œc are the anodic and cathodic transfer coefficients, Cyer,; are 
the reference reactant molar concentrations, a is the volumet- 
ric specific surface area of catalyst (per unit volume of catalyst 
layer), nj, j € {a,c}, is the overpotential and R is the universal 
gas constant (8.314 Jmol~! K~!). The quantity a is a func- 
tion of the mass specific Pt surface area (Pt surface area per 
unit mass of Pt), apt, the Pt loading, me, and the catalyst-layer 
thickness, L 


a ptMpt 
am E. 
L 


The overpotentials, nj, are defined precisely through the rela- 
tionships 


Ne = Y — ġ — Eo, Na = Y — 9, (39) 
where Eo is the so-called open circuit potential, assumed con- 
stant. The value of reactant concentrations in the Butler -Molmer 
expression, (38), ought to reflect accurately that reaction takes 
place at the surfaces of the carbon agglomerates. The concentra- 
tions at these surfaces, Cie i € {O, H}, are generally different 
from the bulk values in the electrolyte, particularly given the 
restricted diffusion of reactants. The Cf, can be related to the 
bulk values, C; e, by balancing the rate of reaction with the rate 
of diffusion of reactant to the surface of the agglomerates (at 
steady state). This mass balance can be approximated as follows 


s Vj s 
V(Cie — Cie) = Kaz, T, Cie), (40) 


where y is the rate of oxygen diffusion through the elec- 
trolyte/water film to the surface of the agglomerates (in s7!). 
Using the definition of R ;(n j, T, C} e) and solving the resulting 
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Table 5 

Heat, charge and mass transfer/transport properties in the base case 

Symbol Quantity Size References 
Dp,o Oxygen diffusivity in free space 3.8 x 107° (e(1 — s)T)?’/? Pm? s7! [35] 
Dan Hydrogen diffusivity in free space 107° (e(1 — s)T)?/?/ P m? s7! [35] 
Dy Vapour diffusivity in free space 4.1 x 1079 (e1 — s)T 3/2 /P m? s7! [35] 
Dy Nitrogen diffusivity in free space 3.6 x 107° (e(1 — s)T)/? /P m? s7! [35] 
Dı Oxygen diffusivity in liquid water (60°C) 4.82 x 107°? m? s7! [35] 
De H Hydrogen diffusivity in nafion 1 x 107!0 m? s7! 

hpe,j Oxygen mass transfer rate 10° s7! est. 
Ho Hy O2(H2) Henry’s law constant (dimensionless) 0.3 (0.6) [53] 
Ka (Kd) Absorption (desorption) constant 1075/9 (1075/3) ms) [2] 
Kj (kg) Absolute permeability of CCL (GDL) 107! (8.7 x 107!) m? [54,47] 
Mm Liquid water viscosity 107? kgm! s7! [4] 
a Surface tension 0.07Nm7! [4] 
k Thermal conductivity of catalyst layers 0.67 W m~! K7! [51] 
km Thermal conductivity of membrane 0.67 W m7! K7! [51] 
ke Thermal conductivity of GDL 1.67 W m~! K7! [51] 
PC Heat capacitance of water 4.187 x 106 J m~? K7! 

PgCg Heat capacitance of air 10? Jm-3 K7! 

PmCm Heat capacitance of membrane 2.18 x 10 J m~’, KT! 

PeCe Heat capacitance of carbon 1.61 x 106 J m~? K7! 

— Ac Entropy associated with ORR 163.7J mol! K7! [52] 
— ^a Entropy associated with HOR? 0Jmol~! K7! [52] 
Os Electronic conductivity before correction 500S m7! [38] 
kdes HO desorption coefficient 100s~! 

Kads HO adsorption coefficient 10s7! 


equation for C}, then yields 


S = IC where 
LEST y + Eer (e%aFn;/RT = e ee fn j/RT)’ 


_ Vjairef,i 


= l (41) 
n FC et? 
so that the final form of the reaction rate is 
Kan, T, CE.) = Rij, T, Cie) 
= n YEerCie (e%aFn;/RT SC eee Fnj/RT) (42) 


vr Y + er (eu! — e~e Fn j/RT) ` 


In order to relate the diffusion rate y to the microscopic 
properties of the catalyst layers, we define an electrolyte film 
thickness, ôe, an agglomerate radius, Rag, and the number of 
agglomerates per unit volume, N 


s Dei s ADg,i 
V(Cie — Cie) ~ AX — (Cie — Ci e), sothat y= ; 
i R be ` i ôe 
KE 
molar flux 
(43) 


where A is the specific surface area of agglomerates (we assume 
that all of the surface area is covered, whereupon A = der RN ). 
The volume of electrolyte attached to the surface of each agglom- 
erate is e£ /N, which, from our assumptions, covers the entire 
surface of the agglomerate. The thickness 45 (without swelling) 
and e are therefore related as follows 


1/3 
SCH R (44) 


The electrolyte swelling results in a volume change equal to (per 
agglomerate) (ef — e$) /N. The electrolyte film thickness is then 
given by 


1/3 


af — e£ 
Je = €0) = Rag. (45) 


GES 

de = | (Rag + doy + AGN 
When s > sx, as previously defined, the liquid water present 
in the catalyst layers will provide additional resistance to the 
oxygen, with a different diffusion coefficient Dı. Assuming that 
the electrolyte is hydrophilic, any liquid water is taken to coat the 
entire surface of the agglomerates. Based on these assumptions, 
the thickness of the water layer is given by 


1/3 
gy = (Ry +85) a) ag + Be) (46) 
LS ag e Aach ag cl, 


with a total film thickness (electrolyte and water) 
ô = be + 5 (47) 


The quantity y then has to be modified; it is decomposed into 
a term arising from diffusion through the water layer, y}, and a 
term of the type described above, ye. To approximate the concen- 
tration of reactant at the agglomerate surfaces two balances need 
to be performed, one for diffusion through the water, to relate 
the bulk concentration to concentration at the water/electrolyte 
interface, and one between diffusion through the electrolyte and 
reaction, to relate the latter concentration to C H e- The reaction 
rate then takes the form (42) with 


YI Ve A DN 


= ———.,_ where = —, = ; 
Vit Ye g A n 


Al = 47 (R ag + 8e) N. (48) 


Y 
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The Henry constant H; would likewise require modification. For 
oxygen, it is approximately 0.3 for water under typical operating 
conditions, and was assigned the value 0.15 for Nafion in [16]. 
However, given the uncertainty in its value we assume that the 
water value holds throughout, both for oxygen and hydrogen. 
Note that the quantity N is a function of both the radius of 
the agglomerates and their distribution, for example, the more 
densely packed the agglomerates the larger N. 
The current flow is the H3O7 flow, so that the current density, 
I' (in A m°), is given by integrating the volumetric value across 
the CCL 
I =4Fr YEeCi,e (ENT IT — eae PnjRT) 
wo Vteer (e%aFn;/RT = e=% Fn;/RT) 


dx. (49) 


2 


2.7. Initial and boundary conditions 


For the discussion of the boundary conditions we recall Fig. 
1. At the interface between the membrane and catalyst layers, 
x = x3 and x = x4, the gas phase fluxes are taken to be zero; that 
is, the gas species are assumed not to penetrate the pore space in 
the membrane. Similarly, we assume that the fluxes of reactants 
in the ionomer phase at the interfaces between the catalyst and 
gas-diffusion (or micro-porous) layers, x = x2 and x = xs, are 
small, i.e., that the reactants enter predominantly through the gas 
phase and that mass exchange becomes significant in the catalyst 
layers. Since protons and dissolved species (reactants and water) 
have no effective transport mechanism in the gas-diffusion (or 
micro-porous) layers, their fluxes at these interfaces are taken to 
be zero. 
Cie db Ke 5 KE 


EE SE 


X = X2, X5 
(50a) 


aCi 
dr 
At the interfaces between the channels and GDL, the gas con- 
centrations and temperature are prescribed, with oxygen only at 
the cathode side and hydrogen only at the anode side. The cell 
voltage, Uc, is measured at the cathode channel/GDL interface, 
and at the anode channel/GDL interface we assume a zero con- 
centration of electrons. 


Ci(0) = Cie (i € {O, N, v}), 
yO) = Uc, wwx7) = 0, 


The values of C ve and C v,a are calculated from the water activity 
in the channels, ge e and aw,a (for cathode and anode respec- 
tively). For the cathode, the activity is defined as Xy cPc / Psat, 
where X v,c is the molar fraction of vapour in the cathode channel, 
P, is the cathode-channel gas pressure and BA is the cathode- 
channel saturation pressure. The saturation pressure, a function 
of temperature, is given by the formula in [37] 


logio Psat = —2.1794 + 0.02953(T — 273) — 9.1837 


=0, ie {O0,H,N, v}, (50b) 


X = X3, X4 


Ci(x7) = Cia (i € {H, N, vii, 
TO) = Te., T(x) = Ta. (50c) 


x10~°(T — 273) + 1.4454 x 107" (T — 273%. 
(50d) 


From these relationships we derive C ve in the cathode as follows 


` RT Paste = 
= e (50e) 
Pe 


where Ĉe is the total gas molar concentration in the cathode 
channel, calculated from the cathode-channel pressure and tem- 
perature 


P. 
E 
RT; 


(S0f) 


The same procedure is applied at the anode channel. 

The final boundary condition is that for the liquid water at the 
interfaces between the gas channels and the GDL. It is common 
in modelling studies, for example [7,6,5,3], to assume either a 
zero saturation or zero liquid water flux at this location, or along 
portions of the channel/GDL interface in two dimensions. In the 
approach of Weber and Newman [43] the liquid water flux is 
assumed zero if the gas pressure is greater than the liquid pres- 
sure, otherwise the capillary pressure is assumed zero. However, 
as with the zero flux conditions, this implicitly assumes that the 
strength of the flow in the channel has no impact on the liquid 
water levels at the GDL surface. In [44] Meng and Wang specify 
the saturation at the GDL/channel interface, which, as they cor- 
rectly state, depends sensitively on the gas flow in the channel, 
current density and wettability [45]. 

The accumulation and removal of liquid water along the chan- 
nel is a complicated process. Water is expelled from the GDL 
through preferential openings and forms droplets attached to the 
surface. These droplets can grow or coalesce to a form larger 
droplets comparable in size to the channel dimensions, which 
results is the formation of a liquid film [46]. Ultimately, the film 
is wicked along the channel walls toward the exit. The greater 
the flow velocity in the channel, the more effective the removal 
of liquid water. A detailed model of this process, accounting for 
the surface properties of the GDL, phase change, surface ten- 
sion, two dimensional gas flow and hydrophillicity of the channel 
walls, is not currently possible. To approximate the physics at 
the interfaces we adopt the following steady-state flux condition 
at x = 0 and x = x7 


Sige, Je 
ax p=y, J a,Cy, 


(50g) 
where &2; = 0 corresponds to zero water removal from the anode 
channel, j = a, or cathode channel, j = c. This form is moti- 
vated in two ways. Firstly, there should be no removal at zero 
saturation. Secondly, the parameter (2; can be seen naturally as 
(in an average sense) the reciprocal of a water film thickness on 
the surface of the GDL. When this thickness approaches zero, the 
saturation is physically zero, which is expressed by 1/Qj; — 0 
in Eq. (50g). The thickness, and therefore Q ;, is likely to depend 
sensitively on the flow rate in the channel [46]. 

Initial conditions for pressure, temperature and vapour con- 
centrations will typically be given by the ambient conditions 
in the channels. The water content of the membrane/ionomer is 
assumed to be given by equilibrium with vapour in the channels. 
The saturation at startup is assumed zero. 
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Fig. 2. The left-hand figure shows polarization curves generated from a sweep cycle at 20 mV s~! using the parameter values in Tables 4 and 5, with no micro-porous 
layer. See Fig. 3 for profiles of saturation and membrane water content during the sweep cycle. In the right-hand figure the effect of the sweep rate is illustrated, with 


arrows indicating the locations of the thresholds. 


2.8. Numerical details 


The governing equations and boundary conditions laid out 
above were solved using the finite-element method implemented 
in COMSOL 3.2a. The discretization of the equations and 
boundary conditions was achieved on a uniform grid using 
quartic Lagrange polynomials as trial and test functions, allow- 
ing the number of grid points to be kept small (typically 64 
or 128). The switch functions were substituted with hyper- 
bolic tanh functions to smooth the discontinuities, a standard 
procedure. 

The default set of parameter values used for the base case is 
given in Tables 4 and 5. No fitting parameters were used since we 
are interested in a qualitative picture of the behaviour. It would be 
possible through a single fitting of one or more parameters, such 
as exchange current density or platinum specific surface area 
(whose values are in any case highly uncertain), to obtained a 
closer quantitative match, but there is no loss of information in 
not doing so. 
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3. Results and discussion 
3.1. Validation and sweep rate effect 


To validate the model we make use of the results in [24] 
of ramp sweeps experiments in potentiostatic and galvanostatic 
mode. We attempt to capture the observed trends with specific 
examples of the effects on performance of the channel condi- 
tions and sweep rate. We then discuss results relating to certain 
microstructural properties of the gas-diffusion and micro-porous 
layers. 

The very small active MEA area (1 cm?) and high flow rates 
(200-400 mL min~!) in the experimental study in [24] ensured 
that the cell behaviour was very close to one-dimensional. The 
cell was conditioned with humidified hydrogen and air before 
switching to dry gases. In our results we shall consider a contin- 
uous input of partly humidified gases; simulation of the precise 
conditions is possible but the details of the humidification pro- 
cess are not provided in [24]. Moreover, a qualitative match 
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Fig. 3. Profiles of saturation in the cathode and membrane water content at selected moments during the forward and backward sweeps corresponding to Fig. 2 with 


a sweep rate of 20mVs~!. The corresponding voltages are indicated. 
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Fig. 4. Profiles of saturation in the cathode and membrane water content at selected moments during the forward and backward sweeps corresponding to the example 


with 10mVs~! in Fig. 2. 


is achieved even with the continuously humidified input. Note 
finally that in [22], in which the same experimental procedure 
is presented and modelled, the authors fail to capture the most 
striking feature of a crossing of the polarization curves generated 
from the backward and forward sweeps in potentiostatic mode 
(to be discussed next). This phenomenon is illustrated in Fig. 1 
of [24] and demonstrated in the many examples in both papers. 

We refer first to Fig. 2, which demonstrates the polarization 
curve for the base-case parameters in Tables 4 and 5, without 
a micro-porous layer. What is immediately striking in this plot 
is the intersection of the curves at approximately 0.75 A cm? 
(0.55 V), marking the boundary between worse performance on 
the backward sweep for cell voltage below the value at the cross- 
ing point, and better performance for cell voltage above this 
value—compare this figure with the experimentally generated 
potentiostatic curves in [24,22]. 

To explain the trend we now refer to Fig. 3, showing the evolu- 
tion of the membrane/ionomer water content and the saturation 
in the cathode at selected times during the sweep cycle. The 
corresponding voltages at these times are indicated. It is clear 
that during the backward sweep the saturation level is greater at 
each cell voltage value, a likely consequence of the different time 


—— forward a =0.6 
— - -backward a =0.6 


om forward a =0.3 
w 


e, backward a ,=0-3 


cell voltage V 


"o 025 05 075 1 1.25 
current density A om" 


scales associated with evaporation and condensation. The mass 
transport limitations arising from the presence of liquid water 
are therefore more strongly felt at low cell voltage, as indicated 
in Fig. 2. However, once the mass transport limitations are over- 
come (the polarization curve enters the so-called ohmic regime) 
the greater saturation level of the backward sweep lowers the 
membrane resistance, in relation to the forward sweep. 

Fig. 2 further shows polarization curves for the base case 
parameters with sweep rates of 10 and 40 mV s~!, demonstrating 
that an increase in sweep rate moves the threshold cell voltage 
to a lower value and threshold current density to a higher value. 
These results are qualitatively very similar to Figs. 3 and 4 of 
[24], with a slight discrepancy in comparison to the former fig- 
ure; the thresholds for 5-10 mV s~! in Fig. 3 of [24] are very 
close, and seem to be at odds with Fig. 4 of [24]. 

The trend exhibited in the right-hand plots of Fig. 2 can be 
explained in relation to the response time of the system to water 
production and condensation. The slower the sweep rate the 
longer the time available to produce water in the catalyst layer, 
which leads to increasing saturation levels at 0.1 V (Fig. 4 shows 
the profiles of saturation and membrane water content during the 
sweep cycle at 10mVs~!). The higher saturation levels dictate 
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Fig. 5. The effects on the threshold value of variations in the water activity and temperature in the channels. The arrows indicate the locations of the thresholds. 
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Fig. 6. The effects of variations in the pore size (dg) and the contact angle (68) in the GDL. 


that the time required to evaporate the water is longer, therefore, 
in general, pushing the threshold current density to lower values 
and threshold cell voltage to higher values as the sweep rate 
is decreased. This, however, will be offset in some measure by 
the increased humidification of the membrane on the backward 
sweep when the saturations levels are greater, so that it is quite 
possible for the two processes to cancel each other. 


3.2. Humidification and temperature effects 


The effects of variations in the temperature and activity of 
the channel streams can be seen in Fig. 5. The left-hand figure 
demonstrates that the threshold phenomenon is much more pro- 
nounced for low water activity and in fact almost disappears as 
activity approaches unity (fully humidified). Although details of 
the experiments with humidified streams are not given in [24], 
the authors do state that in their experiments the threshold phe- 
nomenon was not apparent under these conditions. That this 
trend is observed is not surprising since the greater flooding of 
the CCL at higher activity demands that the time to evaporate 
enough liquid water on the backward sweep to transition from 
the mass-transport regime to the ohmic regime is longer, and in 
fact occurs towards the kinetic part of the curve. In this part of 
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the curve the current density is low and the cell voltage is high 
and their values are not dominated by the membrane resistance. 

The right-hand plots in Fig. 5 show the effect of channel 
temperature: raising the temperature in the channels reduces 
the threshold voltage and increases the threshold current den- 
sity. This result is in complete agreement with the potentiostatic 
sweep curves in [24]. The higher temperature in the channels 
causes reduced saturation levels because of a slower condensa- 
tion rate. This results in a greater limiting current density at the 
higher temperature, despite the slight reduction in the reactant 
concentrations in the channels (by the ideal gas law). Conse- 
quently, a smaller volume of water is required to be evaporated 
on the backward sweep to reach the ohmic regime. 


3.3. Effects of pore size, contact angle and micro-porous 
layer 


We now turn or attention to the effects of two microstructural 
properties in the catalyst and gas-diffusion layers. Fig. 6 demon- 
strates the change in performance as the average pore size and 
contact angle in the GDL are varied. 

As the pore size is decreased the saturation level decreases, a 
consequence of the greater capillary pressure gradients (the pore 
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Fig. 7. A comparison of the profiles of saturation in the cathode during the forward sweeps of the two cases in Fig. 6 with GË = 180° and 6È = 110°. 
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Fig. 8. The effects of a MPL, the right-hand figure showing saturation profiles with and without an MPL using the values in Tables 4 and 5. 


size determines, at fixed porosity, the permeability according 
to the Kozeny—Carman law (5)). Thus the current density at 
0.1 V is greater. The threshold point for dg = 1 pm occurs at 
a higher current density because of the decreased saturation, 
though the difference is small because of the competing process 
of membrane humidification. 

A decreasing contact angle in the GDL leads to an increasing 
saturation because of its decreasing hydrophobicity (see Fig. 7). 
The right-hand plots in Fig. 6 show that in such cases the higher 
saturation forces the threshold point to lower current density 
and higher cell voltage. A competing effect, compounding the 
effect of increased membrane humidification, will be the greater 
retention of water on the backward sweep at 65 = 110°. 

Finally, we review one result relating to the effect of a MPL, 
shown in Fig. 8. The latter compares the sweep cycles with and 
without an MPL, where the total thickness of the gas diffusion 
medium (including the MPL) is kept constant at 200 um and the 
MPL length is 30 um. The MPL can have a dramatic effect on 
the saturation levels in the cell. In general, the higher capillary 
pressure gradients in the smaller pores of the MPL lead to lower 
saturation levels, as seen in the right-hand plots of Fig. 8. The 
lower saturation levels force the current density at 0 V to a higher 
value. The current density at the threshold point is only slightly 
greater with the MPL than without, despite the very different 
levels of saturation on the two sweeps. 

As with Fig. 6 it appears that with severe flooding in one 
case the increase in membrane hydration (and therefore con- 
ductivity) on the backward sweep can offset the flooding effect 
on the relative locations of the thresholds. Moreover, significant 
deviations in the threshold values are more easily seen when the 
forward and backward sweeps differ visibly in the ohmic part of 
the curve, as in Figs. 2 and 5. 


4. Conclusions and future direction 


We have demonstrated that a one-dimensional, transient 
model that fully accounts for liquid water and membrane hydra- 
tion can qualitatively capture complex phenomena observed 
in fuel cell experiments. In particular, the hysteresis that is 


often observed in laboratory experiments has been reproduced 
and explained. Multi-dimensional extensions of the model can 
answer questions related to the dimensions and geometry of the 
channels, anisotropic media and the strengths of the channel 
flows. On the other hand, the qualitative picture can at least be 
captured with the present model, which can therefore form the 
basis of an initial tool to screen MEA designs, combined with 
sweep cycling experiments. The latter are a convenient mea- 
sure of performance, taking far less time than a steady-state 
polarization curve and possessing all of the features of the lat- 
ter, including the competition between flooding and membrane 
hydration. 

Furthermore, the transient aspect of the model provides 
scope for the study of other transient phenomena such as 
fast-transient polarization curves, startup and shutdown (with 
additional kinetic modelling) and freeze-start (with additional 
modelling of the thawing process). These are directions that are 
currently being pursued. 
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